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We  study  localized  modes  in  uniform  one-dimensional  chains  of  tightly  packed  and  uniaxially  compressed 
elastic  beads  in  the  presence  of  one  or  two  light-mass  impurities.  For  chains  composed  of  beads  of  the  same 
type,  the  intrinsic  nonlinearity,  which  is  caused  by  the  Hertzian  interaction  of  the  beads,  appears  not  to  support 
localized,  breathing  modes.  Consequently,  the  inclusion  of  light-mass  impurities  is  crucial  for  their  appearance. 

By  analyzing  the  problem's  linear  limit,  we  identify  the  system’s  eigenfrequencies  and  the  linear  defect  modes. 

Using  continuation  techniques,  we  find  the  solutions  that  bifurcate  from  their  linear  counterparts  and  study 
their  linear  stability  in  detail.  We  observe  that  the  nonlinearity  leads  to  a  frequency  dependence  in  the  ampli¬ 
tude  of  the  oscillations,  a  static  mutual  displacement  of  the  parts  of  the  chain  separated  by  a  defect,  and  for 
chains  with  two  defects  that  are  not  in  contact,  it  induces  symmetry-breaking  bifurcations. 

DOI:  10.11 03/PhysRevE. 80.06660 1  PACS  number(s):  05.45.Yv,  43.25. +  y,  45.70.-n,  46.40.Cd 


I.  INTRODUCTION 

We  investigate  nonlinear  localized  modes  that  result  from 
configurational  heterogeneity  in  granular  crystals.  Accord¬ 
ingly,  our  study  entails  a  confluence  of  three  key  research 
themes;  intrinsic  localization  through  nonlinearity,  wave 
propagation  in  granular  chains,  and  extrinsic  localization 
through  disorder. 

Intrinsic  localized  modes  (ILMs),  also  known  as  discrete 
breathers,  have  been  a  central  theme  for  numerous  theoreti¬ 
cal  and  experimental  studies  over  the  past  two  decades  [1-5]. 
The  original  theoretical  proposal  of  ILMs  in  prototypical  set¬ 
tings  such  as  anharmonic  nonlinear  lattices  [6,7]  and  the  rig¬ 
orous  proof  of  their  existence  under  fairly  general  conditions 
[8]  motivated  numerous  studies  of  such  modes  in  a  diverse 
host  of  applications,  including  optical  waveguides  and  pho- 
torefractive  crystals  [9],  the  denaturation  of  the  DNA  double 
strand  [10],  micromechanical  cantilever  arrays  [11],  nanome¬ 
chanical  resonators  [12],  superconducting  Josephson  junc¬ 
tions  [13],  Bose-Einstein  condensation  [14],  electrical  lat¬ 
tices  [15],  and  more.  In  particular,  ILMs  seem  to  play  a 
significant  role  in  biophysics,  nonlinear  optics,  and  in  con¬ 
densed  matter  physics.  As  shown  in  Refs.  [16],  they  are  ex¬ 
pected  to  contribute  not  only  to  the  dynamics  by  also  to  the 
thermodynamics  of  the  relevant  systems. 

One-dimensional  (ID)  granular  crystals,  which  consist  of 
closely  packed  chains  of  elastically  interacting  particles, 
have  drawn  considerable  attention  during  the  past  few  years. 
This  broad  interest  has  arisen  from  the  wealth  of  available 
material  types  or  sizes  and  the  ability  to  tune  their  dynamic 
response  to  encompass  linear,  weakly  nonlinear,  and  strongly 
nonlinear  regimes  [17-20].  Such  flexibility  makes  granular 
crystals  perfect  candidates  for  many  engineering  applica¬ 
tions,  including  shock  and  energy  absorbing  layers  [21-24], 
actuating  devices  [25],  and  sound  scramblers  [26,27].  Be¬ 
cause  of  these  possibilities,  it  is  crucial  to  investigate  the 
effects  of  defects  (imhomogeneities,  beads  with  different 
masses,  etc.),  which  allow  the  observation  of  interesting 


physical  responses  such  as  fragmentation,  anomalous  reflec¬ 
tions,  and  energy  trapping  [21-24,28-34]. 

It  is  well  known  from  solid-state  physics  that  localized 
vibrations  in  (linear)  lattices  can  arise  from  extrinsic  disorder 
that  breaks  the  discrete  translational  invariance  of  a  perfect 
crystal  lattice  [35,36].  Such  “disorder”  is  also  well  known  to 
introduce  interesting  wave-scattering  effects  [37].  This  phe¬ 
nomenology  arises  in  a  host  of  physical  applications,  includ¬ 
ing  superconductors  [38],  electron-phonon  interactions  [39], 
light  propagation  in  dielectric  superlattices  with  embedded 
defect  layers  [40],  defect  modes  in  photonic  crystals  [41], 
and  optical  waveguide  arrays  [42-44]. 

In  the  present  work,  we  aim  to  investigate  the  confluence 
of  the  preceding  research  themes  by  examining  the  interplay 
of  “disorder”  (which  induces  localized  modes)  and  nonlin¬ 
earity  in  granular  crystals.  We  note  in  passing  that  the  inter¬ 
action  of  impurities  with  solitary  waves  or  a  continuous  os¬ 
cillatory  signal  in  unloaded  (i.e.,  without  precompression) 
monoatomic  chains  has  been  investigated  both  numerically 
[28,31]  and  experimentally  [45].  In  these  studies,  localized 
oscillations  arise  from  the  interaction  of  an  impurity  particle, 
whose  mass  is  lighter  than  that  of  the  other  particles  in  an 
otherwise  homogeneous  chain,  with  either  a  solitary  wave  or 
a  continuous  oscillatory  signal.  However,  these  localized  os¬ 
cillations  are  all  transient,  fading  away  as  soon  as  the  wave 
leaves  the  vicinity  of  the  impurity. 

Here,  by  contrast,  we  examine  long-lived  localized 
breathing  oscillations,  which  form  robust  nonlinear  localized 
modes  (NLMs)  induced  by  the  impurities  in  strongly  com¬ 
pressed  granular  chains.  We  demonstrate  that  their  frequency 
depends  not  only  on  experimental  parameters  such  as  the 
precompressive  force  and  the  constitution  (material  and  size) 
of  the  impurity  bead  but  also  on  the  inherent  nonlinearity  of 
the  system  (i.e.,  on  the  amplitude  of  the  oscillations).  We 
provide  a  detailed  bifurcation  and  dynamical  analysis  of  a 
monoatomic  chain  with  a  single  lower-mass  bead  (which  we 
call  an  “impurity”)  and  extend  our  considerations  to  mono- 
atomic  chains  with  a  pair  of  lower-mass  impurities.  We  show 
that  the  wave  dynamics  in  the  case  of  nearest-neighbor  im- 
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purities,  which  arise  in  a  chain  as  consecutive  lower-mass 
beads  without  other  intervening  particles,  differs  substan¬ 
tially  from  that  in  the  case  of  larger  separations  between 
impurities.  We  focus  on  the  case  of  next-nearest-neighbor 
impurities,  revealing  a  rich  bifurcation  structure  in  which 
asymmetric  branches  of  solutions  emerge  through  spontane¬ 
ous  symmetry  breaking.  We  monitor  the  dynamical  manifes¬ 
tation  of  this  bifurcation  (and  associated  instability)  and  ex¬ 
amine  how  it  is  affected  by  a  potential  initial  asymmetry  in 
the  impurity  masses. 

The  remainder  of  our  presentation  is  structured  as  fol¬ 
lows.  We  start  by  discussing  the  Hertzian  model  of  a  homo¬ 
geneous  (no  impurity)  granular  crystal.  We  then  perform  lin¬ 
ear  and  nonlinear  analyses  first  for  monoatomic  chains  with  a 
single  impurity  and  then  for  monoatomic  chains  with  a  pair 
of  impurities.  Finally,  we  summarize  our  findings  and  present 
some  possible  future  directions. 


II.  MONOATOMIC  GRANULAR  CHAIN 


The  interaction  between  two  adjacent  elastic  spheres  is 
well  known  to  be  described  by  Hertz’s  law  [46].  The  relation 
between  the  force  F0  exerted  on  two  identical  spheres  and 
the  distance  S0  between  their  centers  results  from  geometric 
effects  and  is  given  by  the  nonlinear  relation 

F0  =  A$2,  (1) 


where 

E\2R 

A=W^r 


(2) 


R  is  the  radius  of  the  beads,  E  is  the  material’s  elastic 
(Young’s)  modulus,  and  v  is  the  Poisson  ratio  of  the  bead 
material. 

The  dynamics  of  a  ID  chain  composed  of  beads  of  a 
single  type  (i.e.,  a  monoatomic  chain)  is  thus  described  by 
the  following  system  of  coupled  nonlinear  ordinary  differen¬ 
tial  equations  [17]: 

Miij  =  A\S(j  +  m,_i  -  W;]+2  -  A[<50  +  Uj  -  m, ■+!]+-,  (3) 


where  m;  is  the  displacement  of  the  zth  bead  from  its  equilib¬ 
rium  position  in  the  initially  compressed  chain,  i 
£  {2,  •  ■  ■  ,N—  1},  and  M  is  the  bead  mass.  Both  ends  of  the 
chain  are  free  ( u0=ul,uN=uN+l ),  subjected  to  static  compres¬ 
sion  forces  F0  which  are  causing  the  initial  displacement  S0 
between  adjacent  particles.  The  bracket  [,v]+  of  Eq.  (3)  takes 
the  value  s  if  s> 0  and  the  value  0  if  s<0,  which  signifies 
that  adjacent  beads  are  not  in  contact. 

In  contrast  with  Ref.  [45],  which  considered  unloaded 
chains,  we  investigate  strongly  precompressed  chains,  in 
which  F0  takes  large  values.  Considering  small-amplitude 
displacements  in  comparison  with  the  initial  ones  caused  by 
the  precompressive  force,  so  that 


implies  that  one  can  Taylor  expand  the  forces  in  a  power 
series.  Keeping  the  displacement  terms  to  fourth  order  leads 
to  the  approximate  (iiK2-Ki-K4”)  model 


Miij  =  K2(Uj+l  -  2m,  +  M,_i)  +  K3[(m!+1  -  M,)2  -  (lij_ !  -  m,)2] 

+  K4[(ui+1  -  «i)3  +  («;_ l  -  w,)3],  (5) 

where 

K2=3-AS( \12,  K2  =  -3-AS-0m,  K4=^AS?12.  (6) 

The  equations  of  motion  (5)  are  an  example  of  the  celebrated 
Fermi-Pasta-Ulam  (FPU)  model  [47-49].  If  we  restrict  our 
consideration  to  very  small  amplitudes  and  velocities,  then 
we  can  neglect  all  the  nonlinear  terms  from  the  equations  of 
motion  and  keep  only  the  harmonic  ( K2 )  term.  The  spectral 
band  of  the  ensuing  linear  chain  has  an  upper  cutoff  fre¬ 
quency  of  cjin=  iAK2lm,  which  corresponds  to  the  lattice  vi¬ 
bration  in  which  the  neighboring  particles  oscillate  out  of 
phase. 

One  of  the  remarkable  features  of  gradually  introducing 
nonlinearity  is  that  its  interplay  with  discreteness  leads  to  the 
emergence  of  localized  modes  even  in  the  absence  of  any 
inhomogeneity.  Such  ILM  solutions  are  generic  features  of  a 
large  class  of  Hamiltonian  lattices  that  includes  the  FPU 
model  [4].  Indeed,  ILMs  have  been  studied  extensively  in 
monoatomic  FPU  chains  [50]. 

One  of  the  canonical  mechanisms  for  the  generation  of 
such  nonlinear  localized  modes  is  the  modulational  instabil¬ 
ity  (MI)  of  the  plane  wave  at  the  band  edge.  A  detailed  analy¬ 
sis  of  this  instability  (bifurcation)  shows  that  the  MI  of  the 
upper  cutoff  mode  manifests  itself  when  the  following  in¬ 
equality  holds  (see  Sec.  4.3  of  Ref.  [4]  and  references 
therein): 

3 K2K4  -  4 K\  >  0.  (7) 

Using  Eqs.  (6),  one  can  show  that  Eq.  (7)  does  not  hold  in 
our  setting.  This,  in  turn,  indicates  that  small-amplitude 
ILMs  bifurcating  from  the  upper  band  mode  do  not  exist  in 
monoatomic  granular  crystals.  The  existence  of  dark  discrete 
breathers  or  large-amplitude  ILMs  remains  an  interesting 
open  question  for  future  investigation. 

III.  MONOATOMIC  GRANULAR  CHAIN  WITH  AN 
IMPURITY 

A.  Model 

Consider  a  ID  monoatomic  chain  of  beads  that  contains 
an  impurity  at  the  Arth  site.  Suppose  that  the  impurity  (the  kth 
bead)  is  made  of  the  same  material  as  the  other  particles  but 
has  a  different  radius.  In  particular,  we  consider  a  homoge¬ 
neous  host  chain  that  is  composed  of  spherical  beads  of  non¬ 
magnetic,  type  316  stainless  steel  (which  have  elastic  modu¬ 
lus  £=193  GPa,  Poisson  ratio  m=0.3,  and  density  p 
=  8027.17  kg/m3  [51])  of  radius  R  =  4.76  mm  and  a  spheri¬ 
cal  steel  impurity  bead  with  some  other  radius  r.  We  treat  the 
radius  of  the  impurity  as  a  free  parameter,  though  in  most 
cases  we  will  use  ;  =  0.8/?.  We  also  suppose  that  the  granular 
chain  is  compressed  by  an  experimentally  accessible  static 
force  of  Fq= 22.25  N. 

The  presence  of  the  impurity  bead  in  the  chain  gives  rise 
to  a  mass  defect  and  leads,  in  turn,  to  changes  in  the  force 


066601-2 


LOCALIZED  BREATHING  MODES  IN  GRANULAR... 


PHYSICAL  REVIEW  E  80,  066601  (2009) 


constants  that  destroy  the  discrete  translational  symmetry  of 
the  crystal.  Recalling  that  the  precompressive  static  force  F0 
induces  an  initial  displacement  <50  between  neighboring 
spheres  of  the  same  diameter  and  defining  (5,  to  be  the  dis¬ 
placement  that  it  induces  between  neighboring  spheres  of 
different  diameters,  the  equations  of  motion  (3)  at  sites  k 
—  1,  k,  and  k+  1  become 


MUk_i  -A[[4)+  Uk_  2~  -  A2[<S]  +  uk_i  -  ir.]+2, 

miik  =  A2[Sl  +  uk_\  -  -  A2[<?i  +  uk-  «i+i]+2. 


Miik+ 1  -  A2[<?!  +  uk-  uk+  +  uk+l  -  uk+ 2]+  , 


1/2 


2  E 


R 

2E\  - 

\  2  /  \R  +  r 

Aj  =  ——  ,  A2 


Rr 


1/2 


3(1  -  v2) 


3(1  -v2)  ’ 


(8) 


where  m  and  r  are  the  mass  and  the  radius  of  the  impurity 
bead  while  M  and  R  denote  the  mass  and  the  radius  of  each 
of  the  remaining  beads. 


B.  Harmonic  potential  approximation  (linear  analysis) 

To  understand  the  underlying  linear  spectrum  of  the  prob¬ 
lem,  we  linearize  Eqs.  (3)  and  (8)  about  the  equilibrium  state 
in  the  presence  of  precompression.  This  yields 


mUi  =  —  2m, ■  +  m,+1),  i  $  {k—  l,k,k  +  1}, 


Miik_x  -  C^{uk_2-  uk~ i )  -  C2(uk-i  -  uk)-, 


muk  =  C2(uk_ |  -  2 uk  +  uk+l), 


Miik+\  -  C2(uk-  uk+i)  -  C\(uk+\  -  uk+ 2), 


cx  =  \aXa  c2  =  ^a2s\'2. 


(9) 


Seeking  stationary  solutions  with  frequency  u>,  we  substitute 


Uj  =  vjem'  (10) 

for  all  j  into  Eqs.  (9)  and  obtain  the  following  eigenvalue 
problem: 
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The  eigenvalue  problem  (11)  determines  the  spectrum  of 
the  extended  phonon  excitations  and  of  the  localized  defect 
mode  centered  at  the  impurity  site.  In  general,  the  presence 
of  impurity  beads  can  create  two  types  of  vibrational  modes: 

(1)  Resonance  modes ,  when  the  mass  of  the  impurity  bead 
is  larger  than  the  mass  of  the  other  beads  (i.e.,  when 
m>M). 

(2)  Localized  modes,  in  the  opposite  case  (i.e.,  when 
m<M). 

Each  resonance  mode  possesses  a  frequency  within  the 
range  of  frequencies  that  constitute  the  phonon  band  of  the 
homogeneous  host  crystal  and  has  a  vibration  amplitude  that 
is  larger  in  the  vicinity  of  the  impurity  bead.  Each  localized 
mode,  on  the  other  hand,  has  a  frequency  fimp  that  lies  above 
the  band  of  the  normal-mode  frequencies  of  the  homoge¬ 
neous  host  crystal  and,  as  shown  in  Fig.  1,  has  a  vibration 
amplitude  that  is  large  at  the  impurity  site  but  decreases  very 
rapidly  with  increasing  distance.  Reference  [45]  used 
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FIG.  1.  (Color  online)  (a)  The  normal-mode  frequencies  of  the 
crystal  in  the  presence  of  the  single  lighter-mass  defect.  The  pres¬ 
ence  of  the  impurity  bead  leads  to  the  manifestation  of  a  localized 
mode  (see  inset)  with  frequency  above  the  band  of  the  perfect 
(monoatomic)  host  crystal,  (b)  Numerically  obtained  frequency  fjmp 
(solid  curve)  of  the  localized  mode  as  a  function  of  the  radius  ratio 
r/R  compared  with  the  analytical  prediction  (dashed  curve).  As 
r/R— >1,  one  sees  that  20.67  kHz  (the  horizontal 

dashed  line)  and  that  the  deviation  between  the  solid  and  dashed 
curves  becomes  larger. 

multiple-scale  analysis  to  obtain  the  analytical  approxima¬ 
tion 


2tt  m112 


(12) 


for  fimp,  where  Fm  is  the  maximum  force  magnitude  of  the 
solitary  wave  envelope  moves  over  the  defect  site.  The  ex¬ 
perimental  setup  in  [45]  is  quite  different  from  our  setting 
because  it  has  no  precompression.  As  discussed  in  Ref.  [45], 
the  maximum  force  magnitude  of  this  pulse  acts  as  a  “local 
precompression  force”  as  it  travels,  which  makes  the  system 
weakly  nonlinear  locally  and  results  in  localized  oscillations 
at  the  impurity  site.  By  contrast,  in  the  present  setting,  the 
chain  is  strongly  compressed  by  a  static  force,  which  acts 
globally  in  a  nonlinear  fashion  and  allows  localized  modes  to 
be  maintained  indefinitely  without  further  external  excita¬ 
tions  after  they  are  excited  initially  (by,  e.g.,  an  actuator  or 
the  impact  of  a  striker  particle).  This  distinction  between 
local  and  global  forces  is  the  key  difference  that  leads  to  very 
long-lived  localized  oscillations  around  the  impurity  bead. 
Such  oscillations  can  last  arbitrarily  long  in  principle,  but  in 
laboratory  experiments  the  presence  of  dissipative  effects 


will  eventually  result  in  their  attenuation  [52]. 

As  illustrated  in  Fig.  1(b),  we  find  excellent  agreement 
between  the  analytical  expression  of  Eq.  (12)  for  Fm  =  F0,  and 
the  frequency  of  the  localized  mode  obtained  by  the  eigen¬ 
value  system  (11)  up  to  radii  ratios  ^  —  0.6  (see  also  Fig.  5  of 
[45]).  In  particular,  for  the  material  parameters  and  the  pre- 
compressive  force  indicated  above,  an  impurity  bead  of  ra¬ 
dius  r=0.6R  (for  which  ^  —  0.216)  yields  a  localized  mode 
with  frequency  fimp~  31.76  kHz,  whereas  Eq.  (12)  predicts 
fa~  30  kHz.  On  the  other  hand,  for  r=0.8R  (which  yields 
j[~0.512),  the  eigenvalue  system  (11)  gives  fimp 
—  23.28  kHz  and  Eq.  (12)  predicts  fa  —  20.03  kHz.  To  pro¬ 
vide  additional  context,  we  remark  that  the  upper  cutoff  fre¬ 
quency  of  the  precompressed  homogeneous  host  crystal  is 

given  by  20.67  kHz.  It  is  clear  that  the  ana¬ 

lytical  expression  in  Eq.  (12)  is  expected  to  be  a  good  ap¬ 
proximation  only  for  m  sufficiently  smaller  than  M.  Other¬ 
wise,  one  has  to  use  the  numerically  obtained  frequency  fimp. 
We  also  note  that  as  the  radius  of  the  impurity  bead  (which  is 
smaller  than  the  regular  beads)  becomes  larger,  the  differ¬ 
ence  between  the  frequency  fimp  of  the  impurity-induced  lo¬ 
calized  mode  and  the  upper  cutoff  frequency  becomes 
smaller.  Indeed,  fimp  —> /m  as  r / R— >  1  [as  shown  in  Fig.  1(b)], 
and  the  localized  mode  becomes  concomitantly  more  ex¬ 
tended. 


C.  Continuation  and  stability  analysis 

In  the  previous  section,  we  examined  the  linear  response 
of  the  granular  crystal  in  the  presence  of  a  single  defect.  We 
now  study  Eqs.  (3)  and  (8)  directly  to  examine  the  nonlinear 
behavior  of  the  chain.  From  a  physical  perspective,  this  al¬ 
lows  us  to  examine  the  interplay  of  nonlinearity  and  “disor¬ 
der.”  When  a  perfect  (i.e.,  homogeneous)  nonlinear  lattice 
supports  ILMs,  the  presence  of  impurities  can  drastically 
change  the  properties  of  such  localized  modes.  This  can  re¬ 
sult  in  several  interesting  phenomena — including  the  pres¬ 
ence  of  asymmetric  impurity  modes  in  nonlinear  lattices  with 
a  light-mass  defect  [53]  and  the  existence  of  stable,  nonlin¬ 
ear,  heavy-mass  impurity  modes  [54].  As  we  have  men¬ 
tioned,  our  nonlinear  lattice  does  not  support  ILMs  [see  the 
inequality  (7)],  so  the  aforementioned  phenomena  are  not 
expected  to  be  present.  However,  the  linear  localization  at 
the  defect  in  conjunction  with  nonlinearity  can  result  in  the 
presence  of  robust  NLMs. 

More  specifically,  it  is  important  to  consider  whether  the 
nonlinearity  of  the  chain  can  support  the  existence  of  local¬ 
ized  modes  with  frequencies  /  such  that  fm  <f<f,mp-  In  the 
linear  limit,  the  chain  does  not  support  vibrations  with  fre¬ 
quencies  in  this  regime.  To  answer  this  question  more  gen¬ 
erally,  we  perform  a  parameter  continuation  in  which  we 
start  from  the  linear  localized  mode  and  systematically 
change  (in  small  steps)  the  frequency  from  fimp  toward  the 
upper  cutoff  fm.  For  each  of  the  intermediate  frequencies,  we 
identify  NLMs  to  high  precision,  via  a  Newton  method  in 
phase  space,  using  free  boundary  conditions  and  chains  with 
N=19  beads.  In  order  to  identify  the  relevant  branch  of  so¬ 
lutions,  we  use  as  an  initial  guess  the  localized  impurity- 
induced  mode  [see  the  inset  of  Fig.  1(a)],  as  this  was  ob- 
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FIG.  2.  (Color  online)  (a)  Continuation  diagram  of  the  nonlinear  impurity  modes.  Insets:  the  spatial  profiles  and  the  corresponding 
locations  of  the  Floquet  multipliers  \  in  the  complex  plane  of  the  localized  modes  with  frequencies  fl  =  22.9  kHz,  /2  =  21.65  kHz,  and 
/„,  =  20.67  kHz.  (b)  One  period  of  the  spatiotemporal  evolution  of  the  localized  mode  with  frequency/^  22.9  kHz.  (c)  The  maximum  of  the 
absolute  value  of  the  Floquet  multipliers  as  a  function  of  the  frequency  fb  of  the  nonlinear  impurity  mode. 


tained  from  the  linear  eigenvalue  problem  (11).  The 
momenta  of  all  of  the  sites  can  be  fixed  to  zero,  following 
[55],  due  to  the  time  reversibility  of  the  system.  For  details 
of  this  continuation  method,  see  Ref.  [4]  and  references 
therein. 

In  Fig.  2,  we  show  the  results  of  the  continuation,  which 
allowed  us  to  obtain  localized  solutions  for  all  frequencies 
/e  [fm,fimp).  In  insets  of  Fig.  2(a),  we  show  three  examples 
of  these  solutions;  they  have  frequencies  f\  =22.9  kHz,  f2 
=  21.65  kHz,  and  fm= 20.67  kHz.  We  examined  the  stability 
of  these  localized  modes  by  computing  their  Floquet  multi¬ 
pliers  Xj  to  describe  the  behavior  of  trajectories  near  the  pe¬ 
riodic  solution.  We  show  the  locations  in  the  complex  plane 
of  the  Floquet  multipliers  for  the  three  NLM  profiles  in  in¬ 
sets  of  Fig.  2(a).  As  is  well  known,  if  all  eigenvalues  Xj  have 
unit  magnitude,  then  the  localized  periodic  solution  is  lin¬ 
early  stable.  However,  if  |\,|  >  1  for  some  j,  then  a  perturba¬ 
tion  along  the  corresponding  eigenfunction  ej  grows  by  the 
factor  | Ay|  after  one  period.  In  Fig.  2(b),  we  show  one  period 
of  the  spatiotemporal  evolution  of  the  localized  mode  with 
frequency/!.  In  Fig.  2(c),  we  show  the  absolute  value  of  the 
maximal  eigenvalue,  as  this  is  associated  with  the  instability 
growth  rate. 

We  found  that  the  family  of  the  localized  solutions  exhib¬ 
its  an  oscillatory  instability.  In  general,  oscillatory  instabili¬ 
ties  can  arise  due  to  the  collision  of  either  two  Floquet  mul¬ 
tipliers  associated  with  spatially  extended  eigenvectors  or 
one  multiplier  associated  with  a  spatially  extended  eigenvec¬ 
tor  and  another  associated  with  a  spatially  localized  one.  A 
careful  study  of  the  unstable  Floquet  multipliers  and  the  cor¬ 


responding  eigenvectors  reveals  that  the  oscillatory  instabili¬ 
ties  are  caused  by  the  collision  of  extended  modes  belonging 
to  the  two  arcs  of  overlapping  continuous  phonon  spectrum 
of  the  Floquet  matrix  [56].  During  the  continuation  of  the 
solutions,  typically  3-5  quadruplets  of  eigenvalues  abandon 
the  unit  circle  after  the  collision  but  return  to  it  soon  after¬ 
ward  in  parameter  space.  As  discussed  in  Ref.  [56],  the 
strength  of  such  instabilities  should  depend  on  the  system 
size  (i.e.,  on  the  number  of  beads  in  the  chain).  In  particular, 
when  the  size  of  the  system  is  increased,  the  magnitude  of 
such  instabilities  weakens  uniformly;  simultaneously,  the 
number  of  such  instabilities  increases  with  system  size  due  to 
the  increasing  density  of  colliding  eigenvalues.  Eventually, 
these  instabilities  vanish  in  the  limit  of  an  infinitely  large 
system.  We  observe  that  the  deviations  of  the  unstable  eigen¬ 
values  from  the  unit  circle  are  bounded  by  0.02,  and  numeri¬ 
cal  integration  of  the  nonlinear  impurity  modes  up  to  times 
of  100T  (where  T  is  the  period  of  the  corresponding  mode) 
reveals  their  robustness. 

It  is  relevant  to  also  note  that  two  pairs  of  the  Floquet 
multipliers  are  always  located  at  +1  in  the  complex  plane. 
One  of  them,  often  called  a  “phase  mode,”  describes  a  rota¬ 
tion  of  the  breather’s  aggregate  phase.  The  second  one  arises 
from  the  conservation  of  total  mechanical  momentum,  which 
is  an  integral  of  motion  of  the  nonlinear  chains.  As  one  can 
see  in  the  insets  of  Fig.  2(a),  the  spatial  profiles  of  the  NLMs 
have  interesting  structure.  In  particular,  they  are  character¬ 
ized  by  a  kink-shaped  distortion  of  the  chain  that  is  caused 
by  the  asymmetry  in  the  interaction  potential  (see  Sec.  4.1.4 
of  [4]).  This  asymmetry  is  evident  in  the  K2-K3-K4  ap- 
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FIG.  3.  (Color  online)  Spatiotemporal  evolution  of  the  displacements  of  the  beads  with  initial  conditions  (a)  uk(0)  =  Sl/ 10  and  U;  =  0  for 
all  j=kk  and  (b)  uk(0)  =  S{  and  u,=0  for  all  j  +  k.  For  case  (a),  the  frequency  of  the  oscillations  of  the  impurity  site  is  about  23.25  kHz  and 
max(nj,/  <5j)  =0.065,  while  for  case  (b),  one  obtains  about  22.47  kHz  and  max(nA./<51)«0.5.  Top  insets:  the  displacements  of  the  impurity 
bead  (which  lies  at  the  fcth  site).  Bottom  insets:  magnifications  of  the  spatiotemporal  evolution  near  the  impurity  site.  Observe  that  after  an 
initial  (transient)  stage  of  energy  shedding  in  the  form  of  sound  waves,  a  nonlinear  localized  breathing  mode  forms  in  the  neighborhood  of 
the  defect. 


proximation  of  the  model,  and  it  arises  directly  from  the  fact 
that  Kk  A  0.  The  NLM  can  be  thus  viewed  as  a  localized 
vibration  that  induces  this  kind  of  distortion  in  the  granular 
crystal.  As  one  approaches  the  edge  of  the  phonon  band  (i.e., 
as  /fc —>/„,),  the  NLMs  become  more  extended  and  gradually 
approach  their  extended  (plane  wave)  linear  counterparts  at 
the  upper  band  of  the  linear  spectrum.  In  this  limit,  the  kink¬ 
shaped  distortion  and  the  maximum  of  the  absolute  values  of 
the  Floquet  multipliers  also  increase. 

D.  Excitation  of  nonlinear  impurity  modes 

In  the  previous  section,  we  demonstrated  that  NLMs  exist 
in  the  gap  between  the  band  of  phonon  modes  of  the  perfect 
crystal  and  the  localized  impurity  mode.  This  demonstrates 
that  their  frequency  depends  not  only  on  the  precompressive 
force  and  the  material  parameters  of  the  defect  chain  (as  is 
the  case  in  the  linear  limit)  but  also  on  the  amplitude  dis¬ 
placement  of  the  impurity  bead — in  other  words,  on  the 
strength  of  the  nonlinearity. 

In  this  light,  we  showcase  in  the  present  section  what  we 
believe  is  the  easiest  way  to  observe  these  NLMs.  Our 
method  is  based  on  the  use  of  a  simple  localized  initial  ex¬ 
citation.  At  time  t= 0,  we  displace  the  bead  at  impurity  site  k 
by  an  amount  that  is  strong  enough  to  ensure  that  the  non¬ 
linear  terms  are  no  longer  negligible  in  the  corresponding 
equations  of  motion.  Meanwhile,  we  keep  all  of  the  remain¬ 
ing  sites  at  rest.  We  then  integrate  the  equations  of  motion 
(3)  and  (8)  using  a  fourth-order  Runge-Kutta  numerical 
scheme.  We  expect  part  of  the  initially  localized  energy  ex¬ 
citation  to  spread  among  the  other  sites.  In  order  to  avoid 
backscatter  of  emitted  waves,  we  consider  a  chain  with  N 
=  500  beads. 

We  generated  long-lived  NLMs  using  two  proof-of- 
principle  simulations.  In  the  first,  we  considered  a  relatively 
weak  initial  displacement  of  the  impurity  site,  uk(0)=Sl!  10, 
which  we  show  in  Fig.  3(a).  In  the  second,  we  considered  a 
relatively  strong  initial  displacement  of  the  impurity  site. 


M/r.(0)=^1,  which  we  show  in  Fig.  3(b).  In  both  cases,  we 
used  the  precompression  and  material  parameters  discussed 
previously.  Observe  that  for  the  strong  initial  displacement, 
the  nonlinearity  causes  a  substantial  distortion  of  the  chain. 
The  frequency  of  the  oscillations  of  the  impurity  site  was 
about  23.25  kHz  and  max(wjt/ <5j) —  0.065  for  the  first  simu¬ 
lation.  For  the  second  simulation,  we  observed  a  frequency 
of  about  22.47  kHz  and  max(uk/  S{)  —  0.5.  Both  sets  of  re¬ 
sults  are  in  excellent  agreement  with  the  continuation  analy¬ 
sis  discussed  above,  which  gave  frequencies  of  23.26  and 
22.41  kHz  for  these  particular  values  of  ma x(ukl  Sk).  These 
simulations  therefore  clearly  illustrate  the  excitation  of  the 
previously  analyzed  NLMs. 


IV.  MONOATOMIC  GRANULAR  CHAIN  WITH  TWO 
IMPURITIES 

A.  Harmonic  potential  approximation  (linear  analysis) 

We  now  consider  a  crystal  with  two  impurities,  which  are 
located  at  sites  k  and  /.  As  before,  we  begin  with  an  analysis 
of  the  linearized  spectrum.  We  first  examine  the  case  of  two 
identical  impurities,  for  which  we  consider  impurity  beads 
made  of  the  same  material  (stainless  steel)  as  the  ones  of  the 
host  chain  (which  is  again  a  perfect  crystal)  but  with  smaller 
radius  (ri  =  r2  =  r=0.8R).  The  left  panel  of  Fig.  4  shows  the 
frequencies  that  correspond  to  localized  modes,  which  we 
obtained  using  the  harmonic  approximation,  as  a  function  of 
the  distance  between  the  impurities.  Interestingly,  for  this 
value  of  radius  ratio,  when  the  impurity  beads  are  in  contact 
(l-k=  1),  the  phonon  spectrum  has  just  a  single  localized 
mode  with  frequency  flmp~ 25.28  kHz.  As  shown  in  inset 
(a)  in  the  left  panel  of  Fig.  4,  the  corresponding  mode  is 
antisymmetric.  By  examining  the  phonon  spectrum  for  the 
case  l-k=  1  as  a  function  of  the  radius  ratio  r/R,  we  found 
that  for  r/R  <0.8  a  symmetric  mode  leaves  the  phonon  band 
and  becomes  progressively  more  localized  as  the  radius  ratio 
is  decreased.  In  particular,  for  r/R  =  0.4  the  phonon  spectrum 
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FIG.  4.  (Color  online)  Left  panel:  the  frequencies  of  the  localized  modes  generated  by  two  identical  impurities  as  a  function  of  the 
distance  (number  of  sites)  between  them.  The  insets  show  corresponding  localized  modes  for  (a)  l—k=  1,  for  which  the  impurities  are  in 
contact  and  [(b), (c)]  l—k= 2.  As  the  impurities  are  placed  increasingly  far  apart,  we  see  that  /i-/2— >0.  Right  panel:  the  normal-mode 
frequencies  of  a  granular  crystal  with  two  impurities  of  different  radii  (r]=0.775R  and  r2=0.8R)  and  l-k=  2. 


contains  two  localized  modes:  an  antisymmetric  one  with 
fi~ 62.26  kHz  and  a  symmetric  one  with  /2  —  38.48  kHz. 
When  /-£> 2,  there  are  two  localized  modes  even  for  r/R 
=  0.8.  In  particular,  for  l-k= 2,  the  frequency  of  the  symmet¬ 
ric  mode  [inset  (b)  in  the  left  panel  of  Fig.  4]  is 
fl~ 23.98  kHz,  and  the  frequency  of  the  antisymmetric 
mode  [inset  (c)  in  the  left  panel  of  Fig.  4]  is 
/2 —  22.14  kHz.  As  the  two  impurities  are  placed  farther 
apart,  fx  -f2  —> 0  and  fi,f2^fimp~  23.28  kHz,  which  is  the 
frequency  of  a  single-impurity  localized  mode  (see  the  dis¬ 
cussion  in  the  previous  section). 

We  now  consider  the  phonon  spectrum  for  the  case  of  two 
different  impurities  with  separation  l-k= 2  and  bead  radii 
r j  =  0.775/?  and  r2=0.8R.  (As  before,  the  impurity  beads  are 
made  of  the  same  material  as  those  in  the  host  chain.)  As  one 
can  see  by  comparing  the  insets  in  the  left  and  right  panels  of 
Fig.  4,  even  a  slight  difference  in  the  radii  of  impurity  beads 
results  in  an  asymmetric  modification  of  the  corresponding 
localized  modes.  In  this  case,  we  also  obtain  slightly  larger 
mode  frequencies:  /]  — 24.45  kHz  and  /2  — 22.44  kHz. 

B.  Continuation  and  stability  analysis  for  two  identical 
impurities 

We  show  the  results  of  our  parameter  continuation  for  the 
case  of  two  identical  impurity  beads  in  contact  ( l—k=  1)  with 
radius  ratio  rlR  =  0.8  in  Fig.  5.  We  find  essentially  the  same 
phenomenology  as  we  obtained  for  granular  crystals  with  a 


single  impurity.  That  is,  we  obtain  a  family  of  weakly  oscil¬ 
latory,  unstable  localized  solutions  due  to  finite  size  effects. 
The  magnitudes  of  the  deviations  of  the  unstable  eigenvalues 
from  the  unit  circle  are  smaller  than  0.025.  Again  as  before, 
the  modes  become  wider  and  the  characteristic  kink-shaped 
distortion  of  the  chain  becomes  larger  as  one  approaches  the 
frequency  fm.  We  show  three  examples  of  this  family  of  lo¬ 
calized  solutions  (with  frequencies  /]  =24.75  kHz, 
/ 2= 22. 55  kHz,  and  /„,  =  20.65  kHz)  in  the  insets  of  Fig. 
5(a).  In  the  rest  of  the  insets,  we  show  the  locations  in  the 
complex  plane  of  their  corresponding  Floquet  multipliers. 

Now  consider  a  granular  crystal  with  two  impurities  that 
are  not  in  contact.  More  specifically,  we  focus  on  the  proto¬ 
typical  case  of  I-k= 2  and  r/ A  =  0.8.  As  indicated  above,  the 
corresponding  phonon  spectrum  contains  two  localized 
modes:  a  symmetric  one,  shown  in  inset  (b)  of  the  left  panel 
of  Fig.  4,  and  an  antisymmetric  one,  shown  in  inset  (c)  of  the 
left  panel  of  Fig.  4.  In  Fig.  6,  we  show  continuation  dia¬ 
grams,  which  display  the  frequencies  of  NLMs  in  terms  of 
the  maximum  displacement  of  one  of  the  impurity  beads 
(/= 4 1 ,  in  a  chain  with  A'' =79  beads)  normalized  to  the  char¬ 
acteristic  (precompression-induced)  displacement  <5j.  First, 
we  examine  the  family  of  solutions  that  arises  from  the  sym¬ 
metric  linear  mode  at/]  —  23.98  kHz  [see  Fig.  6(a)].  Stabil¬ 
ity  analysis  demonstrates  the  presence  of  a  weak  oscillatory 
instability,  which  we  also  observed  for  a  single  impurity  and 
a  pair  of  in-contact  impurities.  Now  consider  the  NLMs  that 
bifurcate  from  the  antisymmetric  linear  mode,  for  which 


FIG.  5.  (Color  online)  (a)  Continuation  diagram  of  the  nonlinear  impurity  modes  for  the  case  of  two  identical  impurity  beads  in  contact 
{l-k=  1).  One  set  of  insets  shows  the  profiles  of  the  localized  modes  with  frequencies  /;  =  24.75  kHz,  f2= 22.55  kHz,  and 
/,„  =  20.65  kHz,  and  the  other  set  shows  the  corresponding  locations  in  the  complex  plane  of  their  Floquet  multipliers  X.  (b)  The  maximum 
of  the  absolute  values  of  the  Floquet  multipliers  as  a  function  of  the  frequency  fb  of  the  nonlinear  impurity  mode. 
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FIG.  6.  (Color  online)  (a)  Continuation  diagram  for  the  nonlinear  impurity  modes  that  originate  from  the  symmetric  linear  mode  of  a 
granular  crystal  with  two  impurities  that  are  separated  by  one  bead  (i.e.,  l-k= 2).  One  set  of  insets  shows  the  profiles  for  the  localized  modes 
with  frequencies  /=  23.68  kHz,  /=  22.26  kHz,  and  /=  20.66  kHz.  The  other  set  shows  the  locations  in  the  complex  plane  of  their  corre¬ 
sponding  Floquet  multipliers  A.,  (b)  Continuation  diagram  for  the  antisymmetric  solution  branch  (Aj),  showing  the  generation  of  asymmetric 
branches  (A2  and  A3)  that  originate  from  the  bifurcation  at /— 21.44  kHz.  The  insets  show  the  wave  profiles  (and  their  Floquet  multipliers) 
of  the  antisymmetric  localized  mode  before  the  bifurcation  (fb=  21.94  kHz),  the  antisymmetric  mode  after  the  bifurcation 
(fb= 21.11  kHz),  and  the  asymmetric  mode  with  frequency  fb= 21.11  kHz.  (c)  The  maximum  of  the  absolute  values  of  the  Floquet 
multipliers  in  terms  of  the  frequency  of  the  nonlinear  impurity  mode  fb  for  the  symmetric  branch,  (d)  Same  as  panel  (c),  but  for  the 
antisymmetric  (Aj)  and  asymmetric  (A2,A3)  branches. 


/2~22.14  kHz.  This  family  of  solutions,  which  corresponds 
to  the  branch  A1  of  the  continuation  diagram  in  panel  (b)  of 
Fig.  6,  is  initially  weakly  unstable  (due  to  the  finite-size  ef¬ 
fects  discussed  previously).  At /— 21.56  kHz,  a  pair  of  Flo¬ 
quet  multipliers  leaves  the  phonon  bands.  The  corresponding 
eigenmodes  are  symmetric  and  become  progressively  more 
localized  as  the  frequency  decreases  [57].  At /— 21.44  kHz, 
these  two  localized  modes  collide  at  the  (+1,0)  point  on  the 
unit  circle,  giving  rise  to  a  strong  instability  (called  a  “har¬ 
monic  instability”)  that  is  connected  to  a  bifurcation  of  the 
corresponding  NLM.  Two  new  families  of  solutions 
(branches  A2  and  A3)  emerge  from  this  bifurcation.  Aside 
from  the  kink-shaped  distortion  of  the  system,  these  families 
are  symmetric  with  respect  to  each  other  and  weakly  un¬ 
stable  due  to  finite-size  effects.  Thus,  this  bifurcation  is 
somewhat  reminiscent  of  a  pitchfork  bifurcation  [58].  In  the 
case  of  the  newly  formed  branches  past  the  bifurcation  point 
(i.e.,  A 2  and  A3),  and  particularly  at/~21.34  kHz,  the  lo¬ 
calized  pair  of  eigenmodes  that  has  formed  enters  the  band  of 
eigenvalues  associated  with  extended  perturbations,  giving 
rise  to  a  new  oscillatory  instability.  Although  this  kind  of 
instability  is  size-dependent,  it  persists  even  in  the  limit  of  an 
infinite  system  [59]  (in  contrast  to  the  oscillatory  instability 
caused  by  the  collision  of  two  extended  modes). 

It  is  worth  noting  that  the  setting  of  granular  chains  with 
two  next-nearest-neighbor  impurities  (i.e.,  with  l-k= 2)  is 
reminiscent  of  double-well  configurations  in  other  contexts. 
For  example,  both  defocusing  and  focusing  nonlinear 


Schrodinger  (NLS)  equations  with  double-well  potentials  are 
known  to  exhibit  “symmetry  breaking”  bifurcations  like  the 
one  discussed  above  [60].  (The  defocusing  case  is  relevant  to 
the  present  setting.)  Such  bifurcations  have  even  been  ob¬ 
served  experimentally  in  both  optical  [61]  and  atomic  sys¬ 
tems  [62]. 

Examining  the  temporal  dynamics  of  the  unstable  anti¬ 
symmetric  mode  evinces  the  symmetry-breaking  phenom¬ 
enon.  To  trigger  the  relevant  instability,  we  use  the  sum  of 
the  unstable  solution  with  /=  2 1.2 14  kHz  and  the  corre¬ 
sponding  unstable  localized  eigenfunction  of  the  lineariza¬ 
tion  around  the  solution  as  an  initial  condition  in  the  full 
nonlinear  equations  of  motion  (3)  and  (8).  Its  dynamical  evo¬ 
lution,  which  we  show  in  Fig.  7(a),  reveals  the  “symmetry 
breaking”  at  f~0.4  ms.  This  is  followed  by  alternating  os¬ 
cillations  between  the  two  impurity  sites  (corresponding  to 
the  asymmetric  modes  on  A2  and  A3).  As  illustrated  in  Fig. 
7(b),  the  dynamical  evolution  of  the  weak  oscillatory  insta¬ 
bility  of  the  asymmetric  modes  is  somewhat  similar.  As  we 
pointed  out  above,  such  dynamics  is  reminiscent  of  theoret¬ 
ical  [60]  and  experimental  [61]  observations  of  the  instability 
manifestation  in  NLS  equations  with  double-well  potentials. 

C.  Continuation  and  stability  analysis  for  two  distinct 
impurities 

As  we  discussed  above,  we  observe  a  pitchfork-like  bifur¬ 
cation,  indicating  the  emergence  of  two  families  of  asymmet- 
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FIG.  7.  (Color  online)  Spatiotemporal  evolution  of  the  bead  displacements  for  (a)  the  harmonically  unstable  antisymmetric  localized 
mode  with /fc  =  21. 214  kHz  (for  which  max|\|  —  1.437)  and  (b)  for  the  weakly  oscillatory  unstable  asymmetric  localized  mode  with  the  same 
frequency  (for  which  max|\|  =  1.032). 


ric  localized  solutions  if  the  two  impurities  in  the  granular 
chain  are  identical.  Employing  the  aforementioned  analogy 
with  NLS  equations,  this  implies  that  the  observed  dynamics 
corresponds  to  that  obtained  in  a  symmetric  “double  well” 
that  can  be  construed  as  having  been  induced  by  the  identi¬ 
cal,  next-nearest-neighbor  impurities.  A  natural  generaliza¬ 
tion  is  to  consider  how  such  phenomenology  changes  when 
these  impurities  are  distinct.  In  this  case,  the  “double  well” 
becomes  asymmetric,  and  (from  bifurcation  theory)  one  ex¬ 
pects  to  see  something  analogous  to  what  is  sometimes 
called  an  imperfect  pitchfork  bifurcation  [58].  Such  a  sce¬ 
nario,  which  has  been  observed  in  NLS  equations  with  asym¬ 
metric  double  wells  [63],  involves  an  asymmetric  perturba¬ 
tion  of  the  pitchfork  structure,  resulting  in  a  saddle-node 
bifurcation  and  an  isolated  branch  of  solutions. 

To  observe  this  breakdown  and  obtain  the  associated  bi¬ 
furcation  picture,  we  consider  the  case  of  two  distinct  impu¬ 
rities  with  slightly  different  radii  (namely,  rt  =  0.775/?  and 
r2  =  0.8R )  on  the  prototypical  next-nearest-neighbor  case  of 
I-k= 2.  Recall  that  we  showed  the  normal-mode  frequencies 
for  this  configuration  in  the  right  panel  of  Fig.  4.  We  show 
the  continuation  and  the  stability  diagrams  for  the  family  of 


solutions  originating  from  f2~ 22.44  kHz  in  Fig.  8.  (We 
omit  the  continuation  of  the  solutions  that  emerge  from  the 
linear  mode  with  frequency  because  it  resembles  the  one 
shown  in  the  left  panel  of  Fig.  6(a).)  As  suggested  above,  this 
amounts  to  a  saddle -node  bifurcation.  The  branches  of  solu¬ 
tions  that  are  analogous  to  the  Al  and  A3  branches  from  Fig. 
6  collide  at  a  critical  value/— 21.1  kHz  and  disappear.  Lin¬ 
ear  stability  analysis  shows  that  one  colliding  branch  consists 
of  strongly  harmonically  unstable  solutions  and  that  the  other 
consists  of  weakly  oscillatory  unstable  solutions.  The  iso¬ 
lated  branch,  which  manifests  itself  because  we  have  broken 
the  perfect  pitchfork  and  arises  from  the  mode  with  fre¬ 
quency  /2  in  the  linear  limit,  is  only  weakly  unstable  instead 
of  exhibiting  the  strong  instability  we  showed  in  Fig.  6.  Once 
again,  this  is  reminiscent  of  the  NLS  phenomenology  ob¬ 
served  in  Ref.  [63]. 

V.  CONCLUSIONS  AND  FUTURE  CHALLENGES 

In  conclusion,  we  have  investigated  the  formation  of  lo¬ 
calized  modes  due  to  the  interplay  of  nonlinearity  and  disor¬ 
der  (i.e.,  the  presence  of  defects)  in  granular  crystals.  While 


FIG.  8.  (Color  online)  (a)  Continuation  diagram  for  the  nonlinear  impurity  mode  that  originates  from  the  linear  mode  with 
f2~ 22.44  kHz  for  the  case  of  two  distinct  impurities  with  slightly  different  radii  (namely,  r!=0.775R  and  r2  =  0.8R)  on  the  prototypical 
next-nearest-neighbor  case  of  l-k=  2.  This  diagram  also  shows  two  additional  families  of  solutions  that  collide  at/«21.1  kHz.  The  insets 
show  the  spatial  profiles  (and  associated  Floquet  multipliers  X.)  of  the  localized  modes  at/fc  =  22.16  kHz  and  fb  =  20.9  kHz.  In  the  latter  case, 
we  show  an  example  from  each  branch,  (b)  The  maximum  of  the  absolute  values  of  the  Floquet  multipliers  as  a  function  of  the  frequency 
of  the  nonlinear  impurity  mode  fb  for  the  different  branches.  Observe  that  the  isolated  branch  that  stems  from  the  linear  impurity  mode  is 
only  weakly  unstable  in  the  frequency  interval  that  we  have  probed,  whereas  the  two  branches  that  arise  via  the  saddle-node  bifurcation  and 
exist  for/<21.1  kHz  encompass  both  a  strongly  unstable  family  of  solutions  and  a  weakly  unstable  one. 
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previous  research  has  been  dedicated  to  the  computational 
and  experimental  identification  of  such  modes  [18,45],  the 
present  work  identifies  such  modes  in  a  numerically  exact 
(up  to  a  prescribed  tolerance)  form  and  offer  a  systematic 
analysis  of  their  linear  stability.  We  have  argued  that  these 
localized  modes  are  likely  to  be  absent  in  monoatomic 
chains,  and  we  have  illustrated  that  the  inclusion  of  even  a 
single  defect  can  produce  a  linear  defect  mode  whose  inter¬ 
play  with  nonlinearity  generates  a  whole  family  of  localized 
modes.  We  have  demonstrated  that  these  waves  are  robust, 
and  it  should  be  straightforward  to  generate  them  in  experi¬ 
ments  for  a  wide  range  of  initial  conditions.  We  extended  our 
single-defect  investigation  to  multiple-defect  settings,  for 
which  we  examined  the  prototypical  situation  of  crystals 
with  two  defect  sites.  We  identified  a  much  richer  phenom¬ 
enology  in  such  situations,  as  we  observed  families  of  both 
near-symmetric  and  near-antisymmetric  states.  We  subse¬ 
quently  demonstrated  that  the  latter  yields  additional  families 
of  asymmetric  solutions  through  pitchfork-like  bifurcations, 
which  we  analyzed  for  defect  pairs  containing  both  identical 
and  distinct  impurities. 

Although  this  paper  presents  some  of  the  first  steps  to¬ 
ward  a  systematic  understanding  and  classification  of  local¬ 
ized  breathing  modes  in  granular  crystals,  there  are  numer¬ 
ous  interesting  future  directions  that  should  be  pursued.  First, 
one  can  consider  progressively  larger  numbers  of  defects.  By 
analogy  with  NLS  settings  [64],  we  expect  this  to  reveal  not 
only  a  wealth  of  additional  phenomenology  (e.g.,  new  fami¬ 
lies  of  waves  and  more  complicated  bifurcation  structures) 
but  also  implications  for  how  things  progress  toward  the 
infinite-defect  limit.  If  all  defects  occur  only  as  next-nearest- 
neighbors,  we  would  obtain  a  perfect  diatomic  (two  species) 
crystal  in  this  limit.  Such  a  crystal  is  naturally  expected  to 
support  intrinsic  localized  mode  solutions  in  the  band  gap 
between  its  acoustic  and  optical  bands.  Indeed,  our  prelimi¬ 
nary  results  indicate  that  such  gap  breather  solutions  do 


arise.  The  detailed  examination  of  such  modes,  and  their 
higher-dimensional  generalizations  (including  possibly  gap 
vortex  modes  [65]),  are  among  our  future  goals. 

A  different  but  related  problem  is  the  study  of  ILMs  in  the 
equilibrium  configurations  of  statically  and  strongly  com¬ 
pressed  chains,  which  consist  of  rigid  elements  that  are 
coupled  by  torsional  springs  [66].  The  coupled  torsional 
springs  of  the  type  discussed  in  Ref.  [66]  are  appropriately 
modeled  by  the  discrete  sine -Gordon  equation,  whose  time- 
independent  form  leads  to  the  standard  map  analyzed  in  that 
work  in  some  detail.  Because  of  the  local  term  in  the  non¬ 
linearity  (which  is  absent  here),  the  linear  spectrum  is 
bounded,  so  that  wmin  <  u>unear<  wmax.  By  contrast,  in  a 
granular  crystal,  one  instead  has  0  <  wHnear<  u>m.  In  the  dis¬ 
crete  sine-Gordon  case,  the  general  theory  of  bifurcations  of 
such  ILMs  [4]  illustrates  that  such  modes  will  exist  even  in  a 
homogeneous  lattice  (i.e.,  one  without  any  defects);  they  bi¬ 
furcate  from  the  lower  edge  of  the  spectrum  (i.e.,  from  oj|nm) 
and  exist  for  frequencies  a><  oj]mn.  This  feature  is  notably 
absent  in  a  granular  crystal  (where  wmin=0).  Note,  however, 
that  should  a  defect  (e.g.,  via  a  different  mass  or  spring  cou¬ 
pling)  be  inserted  such  that  a  linear  frequency  a>imp  >  <wmax 
arises,  then  a  nonlinear  defect  mode  would  be  expected  to 
emerge  in  the  sine-Gordon  lattice  just  as  in  the  granular  crys¬ 
tal  analyzed  herein.  The  similarities  and  differences  between 
these  systems  would  be  an  interesting  area  for  future  inves¬ 
tigations. 
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